#######################################################################################################################
############ Governmental Responses to Terrorism in Autocracies: Evidence from China ##################################
## Philip B.K. Potter, Associate Professor, Department of Politics, University of Virginia
## Chen Wang, PhD Candidate, Department of Politics, University of Virginia
#######################################################################################################################

## Load packages ######################################################################################################
library(expss)
library(foreign)
library(survival)
library(KMsurv)
library(MASS)
library(stargazer)
library(pastecs)
library(dplyr)
library(tidyr)
library(rmarkdown)
library(knitr)
library(ggplot2)
library(survminer)
library(pastecs)
library(pscl)
library(simPH)
library(ggthemes)
library(margins)
library(lubridate)
library(sandwich)
library(lmtest)
library(gridExtra)
######################################################################################################################
## Load data #########################################################################################################

## Set working directory
setwd()

## Load data
d1<- readRDS("data_governmental_response_terrorism.RDS")

## Check variable names and labels
names(d1)

checkvariable <- data.frame(matrix(nrow = 32, ncol = 2))
colnames(checkvariable) <- c("variable name", "labels")
for (x in 1:32) {
  checkvariable[x,1]<- colnames(d1[x])
  checkvariable[x,2]<- var_lab(d1[,x])
}
View(checkvariable)

######################################################################################################################

## Generate Figure-1 #################################################################################################
dtarget<- group_by(d1, year)
dtarget<- summarise(dtarget, target_civilian=sum(target_civilian),
                    target_government=sum(target_government))

dtarget_ts<- data.frame("year"=1990:2014)
dtarget_ts<- apply_labels(dtarget_ts, year="Year of the event")
dtarget_ts<- left_join(dtarget_ts, dtarget, by="year")
dtarget_ts$target_civilian<- ifelse(is.na(dtarget_ts$target_civilian)==TRUE,0,
                                    dtarget_ts$target_civilian)
dtarget_ts$target_government<- ifelse(is.na(dtarget_ts$target_government)==TRUE,0,
                                    dtarget_ts$target_government)

dfigure1<- data.frame("year"=rep(1990:2014, 2),
                      "attack"=c(dtarget_ts$target_civilian, dtarget_ts$target_government),
                      "target_type"=rep(c("Civilian Target", "Government Target"), each=25))

p1<- ggplot(dfigure1,aes(x=year,y=attack,fill=target_type))+
  geom_bar(stat="identity")+
  scale_fill_manual(values=c('gray78','grey41'))+
  xlab("Year")+scale_x_continuous(breaks = seq(1990, 2014, 2))+
  geom_hline(yintercept=0)+
  ylab("Number of Attacks")+ggtitle("Uyghur-Initiated Terrorist Incidents, 1990 -2014")+
  theme(panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        axis.title.y = element_text(margin = margin(t = 1, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
        text = element_text(size=30),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=1, linetype="solid", colour ="black"),
        legend.position=c(0.2, 0.8),legend.key.size=unit(0.8, "cm"),
        legend.key = element_rect(fill="white"))
p1
######################################################################################################################

## Generate Figure-2 #################################################################################################
dcasualty<- group_by(d1, year)
dcasualty$count<- 1
dcasualty<- summarise(dcasualty,death_total=sum(killed, na.rm=TRUE),
                      injur_total=sum(wounded, na.rm=TRUE),
                      attack=sum(count))
dcasualty$death_per_attack<- dcasualty$death_total/dcasualty$attack
dcasualty$injur_per_attack<- dcasualty$injur_total/dcasualty$attack

dcasualty_ts<- data.frame("year"=1990:2014)
dcasualty_ts<- apply_labels(dcasualty_ts, year="Year of the event")
dcasualty_ts<- left_join(dcasualty_ts, dcasualty, by="year")
dcasualty_ts$death_per_attack<- ifelse(is.na(dcasualty_ts$death_per_attack)==TRUE,0,
                                       dcasualty_ts$death_per_attack)
dcasualty_ts$injur_per_attack<- ifelse(is.na(dcasualty_ts$injur_per_attack)==TRUE,0,
                                       dcasualty_ts$injur_per_attack)
dcasualty_ts$death_total<- ifelse(is.na(dcasualty_ts$death_total)==TRUE,0,
                                  dcasualty_ts$death_total)

dfigure2<- data.frame("year"=rep(1990:2014, 3),
                      "count"=c(dcasualty_ts$death_per_attack, dcasualty_ts$injur_per_attack,
                                 dcasualty_ts$death_total),
                      "casualty_type"=rep(c("Number of Deaths per Attack", "Number of Injured per Attack",
                                          "Total Number of Deaths"), each=25))

p2<- ggplot(dfigure2,aes(x=year,y=count,linetype=casualty_type))+
  geom_line(size = 1.2)+
  scale_linetype_manual(values=c("twodash","solid","dotted"))+
  ylab("Count")+xlab("Year")+scale_x_continuous(breaks = seq(1990, 2014, 2))+ 
  ggtitle("Killed and Injured, 1990-2014")+
  theme(panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        axis.title.y = element_text(margin = margin(t = 1, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
        text = element_text(size=30),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        panel.border = element_blank(),
        legend.position=c(0.4, 0.8),legend.key.size=unit(1.8, "cm"),
        legend.key = element_rect(fill="white"))
p2
######################################################################################################################

## Models 1-7 reported in Table 1 ####################################################################################
## Generate centered IVs
d1$cgrowth<- d1$growth-mean(d1$growth)                   # centered GDP growth
d1$cmcci<- d1$cci-mean(d1$cci)                           # centered Consumer Confidence Index
d1$cliindex<- d1$liindex-mean(d1$liindex)                # centered Li-index
d1$cmajor<- d1$majority_freq-mean(d1$majority_freq)      # centered Majority Frequency
d1$ctwe<- d1$majority_g20-mean(d1$majority_g20)          # centered Majority Frequency among G20 countries
d1$csc<- d1$majority_scmember-mean(d1$majority_scmember) # centered Majority Frequency among Security Council Members


## Run the Cox models with centered IVs
cox1<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+cluster(id), data = d1,method="breslow")
summary(cox1)

cox2<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+
               internet_ratio+sensitive+cluster(id), data = d1,method="breslow")
summary(cox2)

cox3<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+
             internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
summary(cox3)

cox4<- coxph(Surv(duration, delta)~ cliindex+cmajor+cliindex:cmajor+
             internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
summary(cox4)

cox5<- coxph(Surv(duration, delta)~ cmcci+cmajor+cmcci:cmajor+
             internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id),
             data = d1,method="breslow")
summary(cox5)

cox6<- coxph(Surv(duration, delta)~ cgrowth+ctwe+cgrowth:ctwe+
             internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id),
             data = d1,method="breslow")
summary(cox6)

cox7<- coxph(Surv(duration, delta)~ cgrowth+csc+cgrowth:csc+
             internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
summary(cox7)
######################################################################################################################

## Generate Figure-3 #################################################################################################
cox3<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
summary(cox3)


pdata <- data.frame(cgrowth=c(mean(d1$cgrowth,na.rm=TRUE)+sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)+sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-sd(d1$cgrowth,na.rm=TRUE)),
                     cmajor=c(mean(d1$cmajor,na.rm=TRUE)-sd(d1$cmajor,na.rm=TRUE),
                               mean(d1$cmajor,na.rm=TRUE)-sd(d1$cmajor,na.rm=TRUE),
                               mean(d1$cmajor,na.rm=TRUE)+sd(d1$cmajor,na.rm=TRUE),
                               mean(d1$cmajor,na.rm=TRUE)+sd(d1$cmajor,na.rm=TRUE)),
                     internet_ratio=mean(d1$internet_ratio),
                     casualty=mean(d1$casualty,na.rm=TRUE),
                     sensitive=mean(d1$sensitive),
                     urban=mean(d1$urban),
                     target_civilian=mean(d1$target_civilian),
                     bombing=mean(d1$bombing)
                     
) 

fit0<- survival::survfit(cox3, newdata=pdata,se.fit=TRUE, conf.int=0.95,
                         conf.type=c("plain"),censor=FALSE, type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 4),
                      Case=rep(c("High Growth & Majority: Low Frequency",
                                 "Low Growth & Majority: Low Frequency",
                                 "High Growth & Majority: High Frequency",
                                 "Low Growth & Majority: High Frequency"), each=6),
                      surv=c(fit0$surv[,1:4]),
                      upper=c(fit0$upper[,1:4]),
                      lower=c(fit0$lower[,1:4]))

p3<- ggplot(survdata,aes(x=Days,y=surv,linetype=Case))+
  geom_line(size=1.2, color="black")+scale_linetype_manual(values=c("solid","dashed","dotted","twodash"),
                                                           labels = c("High Growth & High Majority Frequency", 
                                                                      "High Growth & Low Majority Frequency", 
                                                                      "Low Growth & High Majority Frequency",
                                                                      "Low Growth & Low Majority Frequency"))+
  xlab("Number of days after the attack")+
  ylab("Probability of non-reporting")+ggtitle("")+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2,fill="grey30")+
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = -6)))+
  labs(caption = "Note: Shaded area is the 95% confidence interval.", face = "italic")+
  theme(panel.background = element_rect(fill = "white"), 
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        text = element_text(size=36),legend.title = element_blank(),
        legend.background = element_rect(fill="white", linetype="solid", colour ="black"),
        legend.text = element_text(size=26),
        legend.position="right",legend.key.size=unit(2, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(hjust = 0, vjust=-0.6,face = "italic"),
        plot.margin = unit(c(2,2,2,2),"cm"))

p3
######################################################################################################################

## Generate Figure-4 #################################################################################################
cox4<- coxph(Surv(duration, delta)~ cliindex+cmajor+cliindex:cmajor+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
summary(cox4)


pdata <- expand.grid(cmajor= c(-9, -4,1,6,11,16,21), 
                     cliindex=c(mean(d1$cmcci)-sd(d1$cliindex),mean(d1$cliindex)+sd(d1$cliindex)),
                     internet_ratio=mean(d1$internet_ratio),
                     casualty=mean(d1$casualty,na.rm=TRUE),
                     sensitive=mean(d1$sensitive),
                     urban=mean(d1$urban),
                     target_civilian=mean(d1$target_civilian),
                     bombing=mean(d1$bombing))
ypred <- predict(cox4, newdata=pdata,se=TRUE,type="lp",reference="sample")
yy <- ypred$fit + outer(ypred$se, c(0, -1.96, 1.96), '*')

survdata<- as.data.frame(yy)
survdata$hr<- exp(survdata$V1)
survdata$lower<- exp(survdata$V2)
survdata$upper<- exp(survdata$V3)
survdata$liindex<- rep(c(-9+mean(d1$majority_freq), -4+mean(d1$majority_freq), 1+mean(d1$majority_freq),
                         6+mean(d1$majority_freq),11+mean(d1$majority_freq),16+mean(d1$majority_freq),
                         21+mean(d1$majority_freq)), 2)
survdata$cat<- rep(c("Li-Index: Low","Li-Index: High"), each=7)


p4<- ggplot(survdata,aes(x=liindex,y=hr,linetype=cat))+facet_grid(.~ cat)+
  geom_line(size=1.2)+scale_linetype_manual(values=c("solid","twodash"))+
  xlab("Majority Frequency (non-centered value)")+
  ylab("Hazard Ratio")+ggtitle("")+
  scale_y_continuous(trans='log2', breaks = c(0.5,1,8,128))+
  geom_hline(yintercept=1, linetype="dotted", color = "grey20",size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = -6)))+
  labs(caption = "Note: Shaded area is the 95% confidence interval.", face = "italic")+
  theme(panel.background = element_rect(fill = "white", colour = "black"), 
        text = element_text(size=36),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        legend.text = element_text(size=26),
        legend.position="right",legend.key.size=unit(2, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(hjust = 0, vjust=-0.6,face = "italic"),
        plot.margin = unit(c(2,2,2,2),"cm"))
p4
######################################################################################################################

## Generate Figure-5 #################################################################################################
cox5<- coxph(Surv(duration, delta)~ cmcci+cmajor+cmcci:cmajor+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id),
             data = d1,method="breslow")
summary(cox5)

## Fix Majority Frequency at high
c1<- c(min(d1$cci), mean(d1$cci)-1*sd(d1$cci), mean(d1$cci)-0.5*sd(d1$cci),mean(d1$cci),
       mean(d1$cci)+1*sd(d1$cci), mean(d1$cci)+1.5*sd(d1$cci))
c2<- c1-mean(d1$cci)
c3<- rep(mean(d1$cmajor)+1*sd(d1$cmajor), each=6)
pdata <- data.frame(cmcci=c2,
                     cmajor=c3,
                     internet_ratio=mean(d1$internet_ratio),
                     casualty=mean(d1$casualty,na.rm=TRUE),
                     sensitive=mean(d1$sensitive),
                     urban=mean(d1$urban),
                     target_civilian=mean(d1$target_civilian),
                     bombing=mean(d1$bombing)) 

fit0<- survival::survfit(cox5, newdata=pdata,se.fit=TRUE,conf.int=.95,
                         conf.type=c("plain"),censor=FALSE, type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 6),
                      Case=rep(c1, each=6),
                      surv=c(fit0$surv[,1:6]),
                      upper=c(fit0$upper[,1:6]),
                      lower=c(fit0$lower[,1:6]))


survdata_days21_high<- filter(survdata,Days==21)

survdata_days21_high$Case<- as.numeric(survdata_days21_high$Case)


## Fix Majority Frequency at low
c1<- c(min(d1$cci), mean(d1$cci)-1*sd(d1$cci), mean(d1$cci)-0.5*sd(d1$cci),mean(d1$cci),
       mean(d1$cci)+1*sd(d1$cci), mean(d1$cci)+1.5*sd(d1$cci))
c2<- c1-mean(d1$cci)
c3<- rep(mean(d1$cmajor)-1*sd(d1$cmajor), each=6)
pdata <- data.frame(cmcci=c2,
                    cmajor=c3,
                    internet_ratio=mean(d1$internet_ratio),
                    casualty=mean(d1$casualty,na.rm=TRUE),
                    sensitive=mean(d1$sensitive),
                    urban=mean(d1$urban),
                    target_civilian=mean(d1$target_civilian),
                    bombing=mean(d1$bombing)) 

fit0<- survival::survfit(cox5, newdata=pdata,se.fit=TRUE,conf.int=.95,
                         conf.type=c("plain"),censor=FALSE, type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 6),
                      Case=rep(c1, each=6),
                      surv=c(fit0$surv[,1:6]),
                      upper=c(fit0$upper[,1:6]),
                      lower=c(fit0$lower[,1:6]))

survdata_days21_low<- filter(survdata,Days==21)

survdata_days21_low$Case<- as.numeric(survdata_days21_low$Case)

##Combine them together
survdata<- rbind(survdata_days21_low,survdata_days21_high)
survdata<- dplyr::rename(survdata, cci=Case)
survdata$case<- rep(c("Low Majority Frequency","High Majority Frequency"), each=6)



p5<- ggplot(survdata,aes(x= cci,y=surv,linetype=case))+
  geom_line(size=1)+scale_linetype_manual(values=c("solid","twodash"))+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  xlab("Consumer Confidence Index (non-centered value)")+
  #scale_x_continuous(breaks=c(97.5, 98.5,99.5,100.5,101.5,102.5,103.5))+
  ylab("Probability of non-reporting on the 21st Day \nafter the Incident")+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 28, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)))+
  labs(caption = "Note: Shaded area is the 95% confidence interval.", face = "italic")+
  theme(panel.background = element_rect(fill = "white", colour = "black"), text = element_text(size=36),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        legend.text = element_text(size=26),
        legend.position="right",legend.key.size=unit(2, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(hjust = 0, vjust=-0.6,face = "italic"),
        plot.margin = unit(c(2,2,2,2),"cm"))
p5
######################################################################################################################

## Models with alternative specifications ############################################################################
## Again, first center the key IVs
d1$cnatural<- (d1$natural_disaster_days)-mean((d1$natural_disaster_days))       # centered Natural Disaster
d1$cuschina<- d1$uschinadis-mean(d1$uschinadis)                                 # centered US-China Distance
d1$cexternal<- d1$external_condition-mean(d1$external_condition, na.rm = TRUE)  # centered external condition
d1$cdomestic<- d1$domestic_condition-mean(d1$domestic_condition, na.rm=TRUE)    # centered domestic condition

cox8<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+idp+success+
             internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
summary(cox8)

cox9<- coxph(Surv(duration, delta)~ cgrowth+cuschina+cgrowth:cuschina+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
               data = d1,method="breslow")
summary(cox9)

cox10<- coxph(Surv(duration, delta)~ cnatural+cmajor+cnatural:cmajor+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
                data = d1,method="breslow")
summary(cox10)


cox11<- coxph(Surv(duration, delta)~ cnatural+cuschina+cnatural:cuschina+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
                data = d1,method="breslow")
summary(cox11)

cox12<- coxph(Surv(duration, delta)~ cdomestic+cexternal+
                cdomestic:cexternal+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
                data = d1[which(d1$year>1994),],method="breslow")
summary(cox12)

cox13<- coxph(Surv(duration, delta)~ cdomestic+cexternal+
                cdomestic:cexternal+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
                data = d1[which(d1$year>2000),],method="breslow")
summary(cox13)
######################################################################################################################

## Generate Figure-6 #################################################################################################

## Plot for Model 8
cox8<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+idp+success+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
summary(cox8)

pdata <- data.frame(cgrowth=c(mean(d1$cgrowth,na.rm=TRUE)+1.2*sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-1.2*sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)+1.2*sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-1.2*sd(d1$cgrowth,na.rm=TRUE)),
                     cmajor=c(mean(d1$cmajor,na.rm=TRUE)-1.2*sd(d1$cmajor,na.rm=TRUE),
                               mean(d1$cmajor,na.rm=TRUE)-1.2*sd(d1$cmajor,na.rm=TRUE),
                               mean(d1$cmajor,na.rm=TRUE)+1.2*sd(d1$cmajor,na.rm=TRUE),
                               mean(d1$cmajor,na.rm=TRUE)+1.2*sd(d1$cmajor,na.rm=TRUE)),
                     success=mean(d1$success),
                     internet_ratio=mean(d1$internet_ratio),
                     casualty=mean(d1$casualty,na.rm=TRUE),
                     sensitive=mean(d1$sensitive),
                     urban=mean(d1$urban),
                     target_civilian=mean(d1$target_civilian),
                     bombing=mean(d1$bombing),
                     idp=mean(d1$idp)
) 


fit0<- survival::survfit(cox8, newdata=pdata,se.fit=TRUE, conf.int=0.95,conf.type=c("plain"),censor=FALSE, type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 4),
                      Case=rep(c("High Growth & Majority: Low Frequency",
                                 "Low Growth & Majority: Low Frequency",
                                 "High Growth & Majority: High Frequency",
                                 "Low Growth & Majority: High Frequency"), each=6),
                      surv=c(fit0$surv[,1:4]),
                      upper=c(fit0$upper[,1:4]),
                      lower=c(fit0$lower[,1:4]))

g8<- ggplot(survdata,aes(x=Days,y=surv,linetype=Case))+
  geom_line(size=1)+scale_linetype_manual(values=c("solid","dashed","dotted","twodash"),
                                          labels = c("High Growth \nHigh Majority Frequency", 
                                                     "High Growth \nLow Majority Frequency", 
                                                     "Low Growth \nHigh Majority Frequency",
                                                     "Low Growth \nLow Majority Frequency"))+
  xlab("Number of days after the attack")+
  ylab("Probability of non-reporting")+ggtitle("Model 8")+
  #geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  theme(axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)))+
  theme(panel.background = element_rect(fill = "white", colour = "black"), 
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        text = element_text(size=20),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=1, linetype="solid", colour ="black"),
        legend.text=element_text(size=18),
        legend.position="right",legend.justification = c("left"),
        legend.key.height=unit(4,"line"),
        legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.margin = margin(0.2,0.2,0.2,0.2,"cm"))

g8

## Plot for Model 9
cox9<- coxph(Surv(duration, delta)~ cgrowth+cuschina+cgrowth:cuschina+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
summary(cox9)

pdata <- data.frame(cgrowth=c(mean(d1$cgrowth,na.rm=TRUE)+sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)+sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-sd(d1$cgrowth,na.rm=TRUE)),
                     cuschina=c(mean(d1$cuschina,na.rm=TRUE)-sd(d1$cuschina,na.rm=TRUE),
                                mean(d1$cuschina,na.rm=TRUE)-sd(d1$cuschina,na.rm=TRUE),
                                mean(d1$cuschina,na.rm=TRUE)+sd(d1$cuschina,na.rm=TRUE),
                                mean(d1$cuschina,na.rm=TRUE)+sd(d1$cuschina,na.rm=TRUE)),
                     internet_ratio=mean(d1$internet_ratio),
                     casualty=mean(d1$casualty,na.rm=TRUE),
                     sensitive=mean(d1$sensitive),
                     urban=mean(d1$urban),
                     target_civilian=mean(d1$target_civilian),
                     bombing=mean(d1$bombing)
                     
) 

fit0<- survival::survfit(cox9, newdata=pdata,se.fit=TRUE, conf.int=0.95,conf.type=c("plain"),censor=FALSE, type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 4),
                      Case=rep(c("High Growth & US-China Distance: Low",
                                 "Low Growth & US-China Distance: Low",
                                 "High Growth & US-China Distance: High",
                                 "Low Growth & US-China Distance: High"), each=6),
                      surv=c(fit0$surv[,1:4]),
                      upper=c(fit0$upper[,1:4]),
                      lower=c(fit0$lower[,1:4]))


g9<- ggplot(survdata,aes(x=Days,y=surv,linetype=Case))+
  geom_line(size=1)+scale_linetype_manual(values=c("dashed","solid","dotted","twodash"),
                                          labels = c("High Growth \nHigh US-China Distance", 
                                                     "High Growth \nLow US-China Distance", 
                                                     "Low Growth \nHigh US-China Distance",
                                                     "Low Growth \nLow US-China Distance"))+
  xlab("Number of days after the attack")+
  ylab("Probability of non-reporting")+ggtitle("Model 9")+
  #geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  theme(axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)))+
  theme(panel.background = element_rect(fill = "white", colour = "black"), 
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        text = element_text(size=20),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.8, linetype="solid", colour ="black"),
        legend.text=element_text(size=18),
        legend.position="right",legend.justification = c("left"),
        legend.key.height=unit(4,"line"),
        legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.margin = margin(0.2,0.2,0.2,0.2,"cm"))

g9

## Plot for Model 10
cox10<- coxph(Surv(duration, delta)~ cnatural+cmajor+cnatural:cmajor+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1,method="breslow")
summary(cox10)


pdata <- data.frame(cnatural=c(mean(d1$cnatural,na.rm=TRUE)+0.5*sd(d1$cnatural,na.rm=TRUE),
                                mean(d1$cnatural,na.rm=TRUE)-0.5*sd(d1$cnatural,na.rm=TRUE),
                                mean(d1$cnatural,na.rm=TRUE)+0.5*sd(d1$cnatural,na.rm=TRUE),
                                mean(d1$cnatural,na.rm=TRUE)-0.5*sd(d1$cnatural,na.rm=TRUE)),
                     cmajor=c(mean(d1$cmajor,na.rm=TRUE)-1*sd(d1$cmajor,na.rm=TRUE),
                               mean(d1$cmajor,na.rm=TRUE)-1*sd(d1$cmajor,na.rm=TRUE),
                               mean(d1$cmajor,na.rm=TRUE)+1*sd(d1$cmajor,na.rm=TRUE),
                               mean(d1$cmajor,na.rm=TRUE)+1*sd(d1$cmajor,na.rm=TRUE)),
                     internet_ratio=mean(d1$internet_ratio),
                     casualty=mean(d1$casualty,na.rm=TRUE),
                     sensitive=mean(d1$sensitive),
                     urban=mean(d1$urban),
                     target_civilian=mean(d1$target_civilian),
                     bombing=mean(d1$bombing),
                     growth=mean(d1$growth)
                     
) 

fit0<- survival::survfit(cox10, newdata=pdata,se.fit=TRUE, conf.int=0.95,conf.type=c("plain"),censor=FALSE, type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 4),
                      Case=rep(c("Natural Disaster High & Majority Frequency: Low",
                                 "Natural Disaster Low & Majority Frequency: Low",
                                 "Natural Disaster High & Majority Frequency: High",
                                 "Natural Disaster Low & Majority Frequency: High"), each=6),
                      surv=c(fit0$surv[,1:4]),
                      upper=c(fit0$upper[,1:4]),
                      lower=c(fit0$lower[,1:4]))


g10<- ggplot(survdata,aes(x=Days,y=surv,linetype=Case))+
  geom_line(size=1)+scale_linetype_manual(values=c("twodash","dotted","solid","dashed"),
                                          labels = c("High Natural Disaster \nHigh Majority Frequency", 
                                                     "High Natural Disaster \nLow Majority Frequency", 
                                                     "Low Natural Disaster \nHigh Majority Frequency",
                                                     "Low Natural Disaster \nLow Majority Frequency"))+
  xlab("Number of days after the attack")+
  ylab("Probability of on-reporting")+ggtitle("Model 10")+
  #geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  labs(caption = "Note: Shaded area is the 95% confidence interval.", face = "italic")+
  theme(axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)))+
  theme(panel.background = element_rect(fill = "white", colour = "black"), 
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        text = element_text(size=20),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.8, linetype="solid", colour ="black"),
        legend.text=element_text(size=18),
        legend.position="right", legend.box = "horizontal",legend.justification = c("left"),
        legend.key.height=unit(4,"line"),
        legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(hjust = 0, vjust=-1,face = "italic"),
        plot.margin = margin(0.2,0.2,0.2,0.2,"cm"))

g10


## Plot for Model 11
cox11<- coxph(Surv(duration, delta)~ cnatural+cuschina+cnatural:cuschina+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1,method="breslow")
summary(cox11)

pdata <- data.frame(cnatural=c(mean(d1$cnatural,na.rm=TRUE)+1*sd(d1$cnatural,na.rm=TRUE),
                                mean(d1$cnatural,na.rm=TRUE)-1*sd(d1$cnatural,na.rm=TRUE),
                                mean(d1$cnatural,na.rm=TRUE)+1*sd(d1$cnatural,na.rm=TRUE),
                                mean(d1$cnatural,na.rm=TRUE)-1*sd(d1$cnatural,na.rm=TRUE)),
                     cuschina=c(mean(d1$cuschina,na.rm=TRUE)-1*sd(d1$cuschina,na.rm=TRUE),
                                mean(d1$cuschina,na.rm=TRUE)-1*sd(d1$cuschina,na.rm=TRUE),
                                mean(d1$cuschina,na.rm=TRUE)+1*sd(d1$cuschina,na.rm=TRUE),
                                mean(d1$cuschina,na.rm=TRUE)+1*sd(d1$cuschina,na.rm=TRUE)),
                     internet_ratio=mean(d1$internet_ratio),
                     casualty=mean(d1$casualty,na.rm=TRUE),
                     sensitive=mean(d1$sensitive),
                     urban=mean(d1$urban),
                     target_civilian=mean(d1$target_civilian),
                     bombing=mean(d1$bombing),
                     growth=mean(d1$growth)
) 

fit0<- survival::survfit(cox11, newdata=pdata,se.fit=TRUE, conf.int=0.95,conf.type=c("plain"),censor=FALSE, type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 4),
                      Case=rep(c("Natural Diaster: High & US-China Distance: Low",
                                 "Natural Diaster: Low &  US-China Distance: Low",
                                 "Natural Diaster: High &  US-China Distance: High",
                                 "Natural Diaster: Low &  US-China Distance: High"), each=6),
                      surv=c(fit0$surv[,1:4]),
                      upper=c(fit0$upper[,1:4]),
                      lower=c(fit0$lower[,1:4]))

g11<- ggplot(survdata,aes(x=Days,y=surv,linetype=Case))+
  geom_line(size=1)+scale_linetype_manual(values=c("twodash","dashed","dotted","solid"),
                                          labels = c("High Natural Disaster \nHigh US-China Distance", 
                                                     "High Natural Disaster \nLow US-China Distance", 
                                                     "Low Natural Disaster \nHigh US-China Distance",
                                                     "Low Natural Disaster \nLow US-China Distance"))+
  xlab("Number of days after the attack")+
  ylab("Probability of non-reporting")+ggtitle("Model 11")+
  #geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  theme(axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)))+
  theme(panel.background = element_rect(fill = "white", colour = "black"), 
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        text = element_text(size=20),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=1, linetype="solid", colour ="black"),
        legend.text=element_text(size=18),
        legend.position="right",legend.justification = c("right"),
        legend.key.height=unit(4,"line"),
        legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(hjust = 0, vjust=-1,face = "italic"),
        plot.margin = margin(0.2,0.2,0.2,0.2,"cm"))

g11

## Plot for Model 12
cox12<- coxph(Surv(duration, delta)~ cdomestic+cexternal+
                cdomestic:cexternal+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
                data = d1[which(d1$year>1994),],method="breslow")
summary(cox12)

pdata <- data.frame(cdomestic=c(mean(d1$cdomestic, na.rm = TRUE)+1.5*sd(d1$cdomestic, na.rm = TRUE),
                             mean(d1$cdomestic, na.rm = TRUE)-1.5*sd(d1$cdomestic, na.rm = TRUE),
                             mean(d1$cdomestic, na.rm = TRUE)+1.5*sd(d1$cdomestic, na.rm = TRUE),
                             mean(d1$cdomestic, na.rm = TRUE)-1.5*sd(d1$cdomestic, na.rm = TRUE)),
                     cexternal=c(mean(d1$cexternal, na.rm = TRUE)-1*sd(d1$cexternal, na.rm = TRUE),
                                  mean(d1$cexternal, na.rm = TRUE)-1*sd(d1$cexternal, na.rm = TRUE),
                                  mean(d1$cexternal, na.rm = TRUE)+1*sd(d1$cexternal, na.rm = TRUE),
                                  mean(d1$cexternal, na.rm = TRUE)+1*sd(d1$cexternal, na.rm = TRUE)),
                     internet_ratio=mean(d1[which(d1$year>1994),]$internet_ratio),
                     casualty=mean(d1[which(d1$year>1994),]$casualty,na.rm=TRUE),
                     sensitive=mean(d1[which(d1$year>1994),]$sensitive),
                     urban=mean(d1[which(d1$year>1994),]$urban),
                     target_civilian=mean(d1[which(d1$year>1994),]$target_civilian),
                     bombing=mean(d1[which(d1$year>1994),]$bombing),
                     growth=mean(d1[which(d1$year>1994),]$growth)
) 


fit0<- survival::survfit(cox12, newdata=pdata,se.fit=TRUE, conf.int=0.95,conf.type=c("plain"),censor=FALSE,type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 4),
                      Case=rep(c("Domestic Cooperative & External Hostile",
                                 "Domestic Hostile & External Hostile",
                                 "Domestic Cooperative & External Cooperative",
                                 "Domestic Hostile & External Cooperative"), each=6),
                      surv=c(fit0$surv[,1:4]),
                      upper=c(fit0$upper[,1:4]),
                      lower=c(fit0$lower[,1:4]))



g12<- ggplot(survdata,aes(x=Days,y=surv,linetype=Case))+
  geom_line(size=1)+scale_linetype_manual(values=c("solid","dashed","dotted","twodash"),
                                          labels=c("Domestic Cooperative \nExternal Cooperative",
                                                   "Domestic Cooperative \nExternal Hostile",
                                                   "Domestic Hostile \nExternal Cooperative",
                                                   "Domestic Hostile \nExternal Hostile"))+
  xlab("Number of days after the attack")+
  ylab("Probability of non-reporting")+ggtitle("Model 12")+
  #geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  theme(axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)))+
  theme(panel.background = element_rect(fill = "white", colour = "black"), 
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        text = element_text(size=20),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.8, linetype="solid", colour ="black"),
        legend.text=element_text(size=18),
        legend.position="right",legend.justification = c("right"),
        legend.key.height=unit(4,"line"),
        legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.margin = margin(0.2,0.2,0.2,0.2,"cm"))

g12

## Plot for Model 13
cox13<- coxph(Surv(duration, delta)~ cdomestic+cexternal+
                cdomestic:cexternal+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1[which(d1$year>2000),],method="breslow")
summary(cox13)

pdata <- data.frame(cdomestic=c(mean(d1[which(d1$year>2000),]$cdomestic, na.rm = TRUE)+1*sd(d1[which(d1$year>2000),]$cdomestic, na.rm = TRUE),
                             mean(d1[which(d1$year>2000),]$cdomestic, na.rm = TRUE)-1*sd(d1[which(d1$year>2000),]$cdomestic, na.rm = TRUE),
                             mean(d1[which(d1$year>2000),]$cdomestic, na.rm = TRUE)+1*sd(d1[which(d1$year>2000),]$cdomestic, na.rm = TRUE),
                             mean(d1[which(d1$year>2000),]$cdomestic, na.rm = TRUE)-1*sd(d1[which(d1$year>2000),]$cdomestic, na.rm = TRUE)),
                     cexternal=c(mean(d1[which(d1$year>2000),]$cexternal, na.rm = TRUE)-2*sd(d1[which(d1$year>2000),]$cexternal, na.rm = TRUE),
                                  mean(d1[which(d1$year>2000),]$cexternal, na.rm = TRUE)-2*sd(d1[which(d1$year>2000),]$cexternal, na.rm = TRUE),
                                  mean(d1[which(d1$year>2000),]$cexternal, na.rm = TRUE)+2*sd(d1[which(d1$year>2000),]$cexternal, na.rm = TRUE),
                                  mean(d1[which(d1$year>2000),]$cexternal, na.rm = TRUE)+2*sd(d1[which(d1$year>2000),]$cexternal, na.rm = TRUE)),
                     internet_ratio=mean(d1[which(d1$year>2000),]$internet_ratio),
                     casualty=mean(d1[which(d1$year>2000),]$casualty,na.rm=TRUE),
                     sensitive=mean(d1[which(d1$year>2000),]$sensitive),
                     urban=mean(d1[which(d1$year>2000),]$urban),
                     target_civilian=mean(d1[which(d1$year>2000),]$target_civilian),
                     bombing=mean(d1[which(d1$year>2000),]$bombing),
                     growth=mean(d1[which(d1$year>2000),]$growth)
)


fit0<- survival::survfit(cox13, newdata=pdata,se.fit=TRUE, conf.int=0.95,conf.type=c("plain"),censor=FALSE,type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 4),
                      Case=rep(c("Domestic Cooperative & External Hostile",
                                 "Domestic Hostile & External Hostile",
                                 "Domestic Cooperative & External Cooperative",
                                 "Domestic Hostile & External Cooperative"), each=6),
                      surv=c(fit0$surv[,1:4]),
                      upper=c(fit0$upper[,1:4]),
                      lower=c(fit0$lower[,1:4]))


g13<- ggplot(survdata,aes(x=Days,y=surv,linetype=Case))+
  geom_line(size=1)+scale_linetype_manual(values=c("solid","dashed","dotted","twodash"),
                                          labels=c("Domestic Cooperative \nExternal Cooperative",
                                                   "Domestic Cooperative \nExternal Hostile",
                                                   "Domestic Hostile \nExternal Cooperative",
                                                   "Domestic Hostile \nExternal Hostile"))+
  xlab("Number of Days after the Attack")+
  ylab("Probability of non-reporting")+ggtitle("Model 13")+
  #geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  labs(caption = "", face = "italic")+
  theme(axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)))+
  theme(panel.background = element_rect(fill = "white", colour = "black"), 
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        text = element_text(size=20),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.8, linetype="solid", colour ="black"),
        legend.text=element_text(size=18),
        legend.position="right",legend.justification = c("right"),
        legend.key.height=unit(4,"line"),
        legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.margin = margin(0.2,0.2,0.2,0.2,"cm"))

g13

grid.arrange(g8, g11,g9,g12,g10,g13, nrow = 3)


